1 Aim of the analysis

Analyze transcript level data from TIER-Seq comparison and do some exploratory data analysis.

2 DESeq2 analysis

Read in data: htseq-count files created with usegalaxy.eu

htseq_CDS <- read.delim("input/Galaxy129-[Column_Join_on_data_115,_data_123,_and_others].tabular", header=TRUE, row.names=1)
htseq_TUs <- read.delim("input/Galaxy158-[Column_Join_on_data_146,_data_147,_and_others].tabular", header=TRUE, row.names=1)
zuordnung <- read.delim("input/20210114_zuordnungIn_20210114_transcript-history.csv", header=TRUE)
row.names(zuordnung) <- names(htseq_CDS)[c(5,1,6,2,4,3,10,7,11,8,12,9)]
names(htseq_CDS) <- zuordnung[names(htseq_CDS),]$name
names(htseq_TUs) <- zuordnung[names(htseq_TUs),]$name
coldata <- read.csv("input/20210125_colData.csv", row.names=1)

Do actual analysis

# create DESeq2 data object
ddsMat_CDS <- DESeqDataSetFromMatrix(countData = htseq_CDS,
                                 colData = coldata[names(htseq_CDS),],
                                 design = ~ condition)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
ddsMat_TUs <- DESeqDataSetFromMatrix(countData = htseq_TUs,
                                 colData = coldata[names(htseq_TUs),],
                                 design = ~ condition)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
# run DESeq
ddsMat_CDS <- DESeq(ddsMat_CDS) 
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
ddsMat_TUs <- DESeq(ddsMat_TUs) 
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
write.csv(data.frame(factor=ddsMat_CDS$sizeFactor), file="output/transcript_sizeFactors.csv")
write.csv(data.frame(counts(ddsMat_CDS, normalized=TRUE)), file="output/Transcript_CDS_normalizedCounts.csv")

2.1 Diagnostic Plots

plotDispEsts(ddsMat_CDS, main="CDS comparison", xlab="Mean of Normalized Counts", ylab="Dispersion")

plotDispEsts(ddsMat_TUs, main="TU comparison", xlab="Mean of Normalized Counts", ylab="Dispersion")

pdf(file="output/DESeq2_Plots/CDS/ddsMat_CDS_DispEsts.pdf", width=4.5, height=4.5)
plotDispEsts(ddsMat_CDS, main="CDS comparison", xlab="Mean of Normalized Counts", ylab="Dispersion")
dev.off()
## png 
##   2
pdf(file="output/DESeq2_Plots/TU/ddsMat_TU_DispEsts.pdf", width=4.5, height=4.5)
plotDispEsts(ddsMat_TUs, main="TU comparison", xlab="Mean of Normalized Counts", ylab="Dispersion")
dev.off()
## png 
##   2
p <- PCA_plot(ddsMat_CDS, "RNA Features")
p 

ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS_PCA.pdf", plot=p, width=9, height=9, units="cm")

p <- PCA_plot(ddsMat_TUs, "TU")
p

ggsave("output/DESeq2_Plots/TU/ddsMat_TU_PCA.pdf", plot=p, width=9, height=9, units="cm")

p <- heatmap_plot(ddsMat_CDS, "RNA Features")
p

ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS_heatMap.pdf", plot=p, width=15, height=12, units="cm")

p <- heatmap_plot(ddsMat_TUs, "TU")
p

ggsave("output/DESeq2_Plots/TU/ddsMat_TU_heatMap.pdf", plot=p, width=15, height=12, units="cm")

2.2 Extract Results

Extract results and use changeAnnotation_DESeq2.py and annotation_locusTags_stand13012021.csv to add annotation to tables.

# extract results
CDS_result_dWT_dIF_1h <- results(ddsMat_CDS, contrast=c("condition", "dWT_1h", "dIF_1h")) # dWT/dIF -> higher in dWT: higher log2FC
write_tsv(rownames_to_column(as.data.frame(CDS_result_dWT_dIF_1h[order(CDS_result_dWT_dIF_1h$padj),])), file="output/DESeq2_resultsTables/results_CDS-1h.tsv")

CDS_result_dWT_dIF_0h <- results(ddsMat_CDS, contrast=c("condition", "dWT_0h", "dIF_0h")) # dWT/dIF -> higher in dWT: higher log2FC
write_tsv(rownames_to_column(as.data.frame(CDS_result_dWT_dIF_0h[order(CDS_result_dWT_dIF_0h$padj),])), file="output/DESeq2_resultsTables/results_CDS-0h.tsv")

CDS_result_dWT <- results(ddsMat_CDS, contrast=c("condition", "dWT_1h", "dWT_0h")) # dWT 1h/0h -> higher in 1h: higher log2FC
write_tsv(rownames_to_column(as.data.frame(CDS_result_dWT[order(CDS_result_dWT$padj),])), file="output/DESeq2_resultsTables/results_CDS_dWT.tsv")

CDS_result_dIF <- results(ddsMat_CDS, contrast=c("condition", "dIF_1h", "dIF_0h")) # dIF 1h/0h -> higher in 1h: higher log2FC
write_tsv(rownames_to_column(as.data.frame(CDS_result_dIF[order(CDS_result_dIF$padj),])), file="output/DESeq2_resultsTables/results_CDS_dIF.tsv")

TUs_result_dWT_dIF_1h <- results(ddsMat_TUs, contrast=c("condition", "dWT_1h", "dIF_1h")) # dWT/dIF -> higher in dWT: higher log2FC
write_tsv(rownames_to_column(as.data.frame(TUs_result_dWT_dIF_1h[order(TUs_result_dWT_dIF_1h$padj),])), file="output/DESeq2_resultsTables/results_TUs-1h.tsv")

TUs_result_dWT_dIF_0h <- results(ddsMat_TUs, contrast=c("condition", "dWT_0h", "dIF_0h")) # dWT/dIF -> higher in dWT: higher log2FC
write_tsv(rownames_to_column(as.data.frame(TUs_result_dWT_dIF_0h[order(TUs_result_dWT_dIF_0h$padj),])), file="output/DESeq2_resultsTables/results_TUs-0h.tsv")

TUs_result_dWT <- results(ddsMat_TUs, contrast=c("condition", "dWT_1h", "dWT_0h")) # dWT 1h/0h -> higher in 1h: higher log2FC
write_tsv(rownames_to_column(as.data.frame(TUs_result_dWT[order(TUs_result_dWT$padj),])), file="output/DESeq2_resultsTables/results_TUs_dWT.tsv")

TUs_result_dIF <- results(ddsMat_TUs, contrast=c("condition", "dIF_1h", "dIF_0h")) # dIF 1h/0h -> higher in 1h: higher log2FC
write_tsv(rownames_to_column(as.data.frame(TUs_result_dIF[order(TUs_result_dIF$padj),])), file="output/DESeq2_resultsTables/results_TUs_dIF.tsv")

Execute in directory “DESeq2_resultsTables”: python changeAnnotation_DESeq2.py results_CDS-1h.tsv results_CDS-1h_annotated.tsv python changeAnnotation_DESeq2.py results_CDS-0h.tsv results_CDS-0h_annotated.tsv python changeAnnotation_DESeq2.py results_CDS_dWT.tsv results_CDS-dWT_annotated.tsv python changeAnnotation_DESeq2.py results_CDS_dIF.tsv results_CDS-dIF_annotated.tsv

python TUs_add_info.py results_TUs-1h.tsv results_TUs-1h_annotated.tsv python TUs_add_info.py results_TUs-0h.tsv results_TUs-0h_annotated.tsv python TUs_add_info.py results_TUs_dWT.tsv results_TUs-dWT_annotated.tsv python TUs_add_info.py results_TUs_dIF.tsv results_TUs-dIF_annotated.tsv

2.3 p-Value plots for different comparisons

pvaluePlot(CDS_result_dWT_dIF_0h, "CDS 0h")

pvaluePlot(CDS_result_dWT_dIF_1h, "CDS 1h")

pvaluePlot(CDS_result_dWT, "CDS dWT")

pvaluePlot(CDS_result_dIF, "CDS dIF")

pvaluePlot(TUs_result_dWT_dIF_0h, "TUs 0h")

pvaluePlot(TUs_result_dWT_dIF_1h, "TUs 1h")

pvaluePlot(TUs_result_dWT, "TUs dWT")

pvaluePlot(TUs_result_dIF, "TUs dIF")

## png 
##   2
## png 
##   2
## png 
##   2
## png 
##   2
## png 
##   2
## png 
##   2
## png 
##   2
## png 
##   2

2.4 Filter for non-NA values

CDS_result_dWT_dIF_0h <- subset(CDS_result_dWT_dIF_0h, !is.na(CDS_result_dWT_dIF_0h$padj))
nrow(CDS_result_dWT_dIF_0h)
## [1] 4318
CDS_result_dWT_dIF_1h <- subset(CDS_result_dWT_dIF_1h, !is.na(CDS_result_dWT_dIF_1h$padj))
nrow(CDS_result_dWT_dIF_1h)
## [1] 5293
CDS_result_dWT <- subset(CDS_result_dWT, !is.na(CDS_result_dWT$padj))
nrow(CDS_result_dWT)
## [1] 4968
CDS_result_dIF <- subset(CDS_result_dIF, !is.na(CDS_result_dIF$padj))
nrow(CDS_result_dIF)
## [1] 5510
TUs_result_dWT_dIF_0h <- subset(TUs_result_dWT_dIF_0h, !is.na(TUs_result_dWT_dIF_0h$padj))
nrow(TUs_result_dWT_dIF_0h)
## [1] 3620
TUs_result_dWT_dIF_1h <- subset(TUs_result_dWT_dIF_1h, !is.na(TUs_result_dWT_dIF_1h$padj))
nrow(TUs_result_dWT_dIF_1h)
## [1] 4009
TUs_result_dWT <- subset(TUs_result_dWT, !is.na(TUs_result_dWT$padj))
nrow(TUs_result_dWT)
## [1] 3853
TUs_result_dIF <- subset(TUs_result_dIF, !is.na(TUs_result_dIF$padj))
nrow(TUs_result_dIF)
## [1] 4009

2.5 Create MA and Volcano plots for different comparisons, count features up-, downregulated

2.5.1 CDS

2.5.1.1 CDS 0h

count_up_down(CDS_result_dWT_dIF_0h, foldchange=0.8, padjusted=0.05)
## [1] "number of features down:  20"
## [1] "number of features up:  94"
p <- volcanoPlot_ggplot(as.data.frame(CDS_result_dWT_dIF_0h), foldchange=0.8, padjusted=0.05, text=TRUE)
p
## Warning: Removed 4204 rows containing missing values (geom_text_repel).
## Warning: ggrepel: 50 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-0h_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
## Warning: Removed 4204 rows containing missing values (geom_text_repel).
## Warning: ggrepel: 90 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
p <- MAplot_ggplot(CDS_result_dWT_dIF_0h, foldchange=0.8, y_axis_label = "Log2 fold-change(rne(WT)/rne(Ts))") + ylim(-5.5,+5.5)
p

ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-0h_MAplot.pdf",plot=p, width=15, height=12, units="cm")

2.5.1.2 CDS 1h

count_up_down(CDS_result_dWT_dIF_1h, foldchange=0.8, padjusted=0.05)
## [1] "number of features down:  433"
## [1] "number of features up:  443"
p <- volcanoPlot_ggplot(as.data.frame(CDS_result_dWT_dIF_1h), foldchange=0.8, padjusted=0.05, text=TRUE)
p
## Warning: Removed 4417 rows containing missing values (geom_text_repel).
## Warning: ggrepel: 817 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-1h_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
## Warning: Removed 4417 rows containing missing values (geom_text_repel).
## Warning: ggrepel: 863 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
p <- MAplot_ggplot(CDS_result_dWT_dIF_1h, foldchange=0.8, y_axis_label = "Log2 fold-change(rne(WT)/rne(Ts))") + ylim(-5.5,+5.5)
p

ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-1h_MAplot.pdf",plot=p, width=15, height=12, units="cm")

2.5.1.3 CDS dWT

count_up_down(CDS_result_dWT, foldchange=0.8, padjusted=0.05)
## [1] "number of features down:  360"
## [1] "number of features up:  302"
p <- volcanoPlot_ggplot(as.data.frame(CDS_result_dWT), foldchange=0.8, padjusted=0.05, text=TRUE)
p
## Warning: Removed 4306 rows containing missing values (geom_text_repel).
## Warning: ggrepel: 618 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-dWT_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
## Warning: Removed 4306 rows containing missing values (geom_text_repel).
## Warning: ggrepel: 645 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
p <- MAplot_ggplot(CDS_result_dWT, foldchange=0.8, y_axis_label = "Log2 fold-change(rne(WT) 1h / 0h)") + ylim(-5.5,+5.5) + scale_colour_manual(values=c("DOWN"="black", "UP"="black", "NO"="#d3d3d3ff"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p
## Warning: Removed 1 rows containing missing values (geom_point).

ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-dWT_MAplot.pdf",plot=p, width=15, height=12, units="cm")
## Warning: Removed 1 rows containing missing values (geom_point).

2.5.1.4 CDS dIF

count_up_down(CDS_result_dIF, foldchange=0.8, padjusted=0.05)
## [1] "number of features down:  636"
## [1] "number of features up:  686"
p <- volcanoPlot_ggplot(as.data.frame(CDS_result_dIF), foldchange=0.8, padjusted=0.05, text=TRUE)
p
## Warning: Removed 4188 rows containing missing values (geom_text_repel).
## Warning: ggrepel: 1277 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-dIF_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
## Warning: Removed 4188 rows containing missing values (geom_text_repel).
## Warning: ggrepel: 1311 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
p <- MAplot_ggplot(CDS_result_dIF, foldchange=0.8, y_axis_label = "Log2 fold-change(Rne-Ts 1h / 0h)") #+ ylim(-5,+5)
p

ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-dIF_MAplot.pdf",plot=p, width=15, height=12, units="cm")

2.5.2 TUs

2.5.2.1 TUs 0h

count_up_down(TUs_result_dWT_dIF_0h, foldchange=0.8, padjusted=0.05)
## [1] "number of features down:  8"
## [1] "number of features up:  78"
p <- volcanoPlot_ggplot(as.data.frame(TUs_result_dWT_dIF_0h), foldchange=0.8, padjusted=0.05)
p

ggsave("output/DESeq2_Plots/TU/ddsMat_TUs-0h_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")

DESeq2::plotMA(TUs_result_dWT_dIF_0h)
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)

pdf(file="output/DESeq2_Plots/TU/ddsMat_TUs-0h_MAplot.pdf", width=4.5, height=4.5)
DESeq2::plotMA(TUs_result_dWT_dIF_0h, xlab="Mean of Normalized Counts", ylab="Log2FC(rne-OE/rne-ts)")
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
dev.off()
## png 
##   2

2.5.2.2 TUs 1h

count_up_down(TUs_result_dWT_dIF_1h, foldchange=0.8, padjusted=0.05)
## [1] "number of features down:  327"
## [1] "number of features up:  309"
p <- volcanoPlot_ggplot(as.data.frame(TUs_result_dWT_dIF_1h), foldchange=0.8, padjusted=0.05)
p

ggsave("output/DESeq2_Plots/TU/ddsMat_TUs-1h_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")

DESeq2::plotMA(TUs_result_dWT_dIF_1h)
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)

pdf(file="output/DESeq2_Plots/TU/ddsMat_TUs-1h_MAplot.pdf", width=4.5, height=4.5)
DESeq2::plotMA(TUs_result_dWT_dIF_1h, xlab="Mean of Normalized Counts", ylab="Log2FC(rne-OE/rne-ts)")
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
dev.off()
## png 
##   2

2.5.2.3 TUs dWT

count_up_down(TUs_result_dWT, foldchange=0.8, padjusted=0.05)
## [1] "number of features down:  268"
## [1] "number of features up:  227"
p <- volcanoPlot_ggplot(as.data.frame(TUs_result_dWT), foldchange=0.8, padjusted=0.05)
p

ggsave("output/DESeq2_Plots/TU/ddsMat_TUs-dWT_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")

DESeq2::plotMA(TUs_result_dWT)
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)

pdf(file="output/DESeq2_Plots/TU/ddsMat_TUs-dWT_MAplot.pdf", width=4.5, height=4.5)
DESeq2::plotMA(TUs_result_dWT, xlab="Mean of Normalized Counts", ylab="Log2FC(rne-OE/rne-ts)")
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
dev.off()
## png 
##   2

2.5.2.4 TUs dIF

count_up_down(TUs_result_dIF, foldchange=0.8, padjusted=0.05)
## [1] "number of features down:  492"
## [1] "number of features up:  524"
p <- volcanoPlot_ggplot(as.data.frame(TUs_result_dIF), foldchange=0.8, padjusted=0.05)
p

ggsave("output/DESeq2_Plots/TU/ddsMat_TUs-dIF_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")

DESeq2::plotMA(TUs_result_dIF)
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)

pdf(file="output/DESeq2_Plots/TU/ddsMat_TUs-dIF_MAplot.pdf", width=4.5, height=4.5)
DESeq2::plotMA(TUs_result_dIF, xlab="Mean of Normalized Counts", ylab="Log2FC(rne-OE/rne-ts)")
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
dev.off()
## png 
##   2

3 Exploratory Data Analysis

3.1 Create table how many RNA features of a certain type are affected by differential expression

features <- rtracklayer::import("input/20210217_syne_onlyUnique_withFeat.gff3")
CDS_features <- subset(features, features$type=="CDS")
TUs <- rtracklayer::import("input/Kopf_4091_TUs_combined.gff3")

types <- rep(c("CDS","5UTR", "3UTR", "tRNA", "rRNA", "ncRNA", "asRNA", "CRISPR", "misc"),2)
updown <- c(rep("up",9), rep("down",9))
df_features <- data.frame(cbind(type=types, updown=updown))

types <- c("CDS","5UTR", "3UTR", "tRNA", "rRNA", "ncRNA", "asRNA", "CRISPR", "misc")
for(t in types){
  t_feat <- subset(features, features$type==t)
  
  subset_CDS <- subset(CDS_result_dWT, row.names(CDS_result_dWT) %in% t_feat$locus_tag)
  
  index_up <- which(df_features$type==t & df_features$updown=="up")
  index_down <- which(df_features$type==t & df_features$updown=="down") 
  
  df_features[index_up,"number_feat_overlap"] <- nrow(subset(subset_CDS, subset_CDS$padj<0.05 & subset_CDS$log2FoldChange>0.8)) # count number of features affected
  df_features[index_down,"number_feat_overlap"] <- nrow(subset(subset_CDS, subset_CDS$padj<0.05 & subset_CDS$log2FoldChange<(-0.8))) # count number of features affected
  
  df_features[df_features$type==t,"number_total"] <- length(t_feat)
} 

df_features
##      type updown number_feat_overlap number_total
## 1     CDS     up                 193         3675
## 2    5UTR     up                  25          979
## 3    3UTR     up                   3           29
## 4    tRNA     up                   0           42
## 5    rRNA     up                   1            6
## 6   ncRNA     up                  28          318
## 7   asRNA     up                  52         1071
## 8  CRISPR     up                   0            3
## 9    misc     up                   0           11
## 10    CDS   down                 292         3675
## 11   5UTR   down                  26          979
## 12   3UTR   down                   2           29
## 13   tRNA   down                   5           42
## 14   rRNA   down                   0            6
## 15  ncRNA   down                  10          318
## 16  asRNA   down                  24         1071
## 17 CRISPR   down                   0            3
## 18   misc   down                   1           11
write.csv(df_features, file="output/RNAfeatures_upDown_dWT_0h1hheat.csv")
# prepare data.frame for barplot
types <- rep(c("CDS","5UTR", "3UTR", "tRNA", "rRNA", "ncRNA", "asRNA", "CRISPR", "misc"),2)
updown <- c(rep("up",9), rep("down",9))
df_features <- data.frame(cbind(type=types, updown=updown))

types <- c("CDS","5UTR", "3UTR", "tRNA", "rRNA", "ncRNA", "asRNA", "CRISPR", "misc")
for(t in types){
  t_feat <- subset(features, features$type==t)
  
  subset_CDS <- subset(CDS_result_dWT_dIF_1h, row.names(CDS_result_dWT_dIF_1h) %in% t_feat$locus_tag)
  
  index_up <- which(df_features$type==t & df_features$updown=="up")
  index_down <- which(df_features$type==t & df_features$updown=="down") 
  
  df_features[index_up,"number_feat_overlap"] <- nrow(subset(subset_CDS, subset_CDS$padj<0.05 & subset_CDS$log2FoldChange>0.8)) # count number of features affected
  df_features[index_down,"number_feat_overlap"] <- nrow(subset(subset_CDS, subset_CDS$padj<0.05 & subset_CDS$log2FoldChange<(-0.8))) # count number of features affected
  
  df_features[df_features$type==t,"number_total"] <- length(t_feat)
} 

df_features
##      type updown number_feat_overlap number_total
## 1     CDS     up                 333         3675
## 2    5UTR     up                  29          979
## 3    3UTR     up                   2           29
## 4    tRNA     up                  17           42
## 5    rRNA     up                   0            6
## 6   ncRNA     up                  28          318
## 7   asRNA     up                  31         1071
## 8  CRISPR     up                   3            3
## 9    misc     up                   0           11
## 10    CDS   down                 226         3675
## 11   5UTR   down                  41          979
## 12   3UTR   down                   1           29
## 13   tRNA   down                   0           42
## 14   rRNA   down                   0            6
## 15  ncRNA   down                  24          318
## 16  asRNA   down                 140         1071
## 17 CRISPR   down                   0            3
## 18   misc   down                   1           11
write.csv(df_features, file="output/RNAfeatures_upDown_dWT_dIF_1hheat.csv")

3.2 Functional Enrichment

go_functional_enrichment(CDS_result_dWT, write=TRUE, path_up="output/enrichment/go_enrichment/dWT_0h_1h_up.csv", path_down="output/enrichment/go_enrichment/dWT_0h_1h_down.csv")
##                    ID       Description GeneRatio BgRatio       pvalue
## GO:0006310 GO:0006310 DNA recombination     7/130 27/2459 0.0003587321
##              p.adjust     qvalue
## GO:0006310 0.03694941 0.03694941
##                                                             geneID Count
## GO:0006310 sll0832/sll6059/sll6109/slr0181/slr5010/slr7005/slr8029     7
##                    ID                              Description GeneRatio
## GO:0008137 GO:0008137 NADH dehydrogenase (ubiquinone) activity    13/240
## GO:0048038 GO:0048038                          quinone binding    13/240
## GO:0042773 GO:0042773 ATP synthesis coupled electron transport     9/240
## GO:0031470 GO:0031470                                     <NA>     8/240
## GO:0015977 GO:0015977                          carbon fixation     8/240
## GO:0015986 GO:0015986   ATP synthesis coupled proton transport     6/240
##            BgRatio       pvalue     p.adjust       qvalue
## GO:0008137 19/2459 8.545747e-10 1.076764e-07 9.805120e-08
## GO:0048038 21/2459 5.353053e-09 3.372423e-07 3.070962e-07
## GO:0042773 14/2459 9.008581e-07 3.783604e-05 3.445387e-05
## GO:0031470 12/2459 2.579196e-06 8.124469e-05 7.398221e-05
## GO:0015977 13/2459 6.143525e-06 1.548168e-04 1.409777e-04
## GO:0015986 10/2459 1.220117e-04 2.562246e-03 2.333207e-03
##                                                                                                             geneID
## GO:0008137 sll0223/sll0520/sll0521/sll0522/sll1732/sll1733/slr0261/slr1279/slr1280/slr1281/slr1291/slr2007/slr2009
## GO:0048038 sll0223/sll0519/sll0520/sll0521/sll0522/sll1732/slr0261/slr1279/slr1280/slr1281/slr1623/slr2007/slr2009
## GO:0042773                                 sll0223/sll0522/sll1732/sll1733/slr1291/slr2006/slr2007/slr2008/slr2009
## GO:0031470                                         sll1028/sll1029/sll1030/sll1031/sll1032/slr0009/slr0012/slr0169
## GO:0015977                                         sll1028/sll1029/sll1030/sll1031/sll1032/slr0009/slr0012/slr0169
## GO:0015986                                                         sll1324/sll1325/sll1326/sll1327/slr1330/ssl2615
##            Count
## GO:0008137    13
## GO:0048038    13
## GO:0042773     9
## GO:0031470     8
## GO:0015977     8
## GO:0015986     6
go_functional_enrichment(CDS_result_dIF, write=TRUE, path_up="output/enrichment/go_enrichment/dIF_0h_1h_up.csv", path_down="output/enrichment/go_enrichment/dIF_0h_1h_down.csv")
##                    ID                           Description GeneRatio BgRatio
## GO:0005887 GO:0005887 integral component of plasma membrane    20/280 75/2476
## GO:0004673 GO:0004673     protein histidine kinase activity    13/280 43/2476
## GO:0000155 GO:0000155   phosphorelay sensor kinase activity    15/280 54/2476
##                  pvalue   p.adjust     qvalue
## GO:0005887 0.0001447987 0.02157500 0.01950971
## GO:0004673 0.0005917603 0.03089794 0.02794020
## GO:0000155 0.0006221061 0.03089794 0.02794020
##                                                                                                                                                                     geneID
## GO:0005887 sll0016/sll0094/sll0477/sll0478/sll0648/sll1124/sll1191/sll1192/sll1404/sll1405/sll1600/slr0210/slr0533/slr1147/slr1249/slr1317/slr1400/slr1647/slr1805/slr2045
## GO:0004673                                                         sll0094/sll0790/sll1124/sll1888/sll1905/slr0210/slr0533/slr1147/slr1212/slr1285/slr1400/slr1731/slr1805
## GO:0000155                                         sll0094/sll0790/sll1124/sll1888/sll1905/slr0210/slr0487/slr0533/slr0829/slr1147/slr1212/slr1285/slr1400/slr1731/slr1805
##            Count
## GO:0005887    20
## GO:0004673    13
## GO:0000155    15
##                    ID                              Description GeneRatio
## GO:0042651 GO:0042651                       thylakoid membrane    36/371
## GO:0008137 GO:0008137 NADH dehydrogenase (ubiquinone) activity    13/371
## GO:0019684 GO:0019684           photosynthesis, light reaction    13/371
## GO:0048038 GO:0048038                          quinone binding    13/371
## GO:0042773 GO:0042773 ATP synthesis coupled electron transport    10/371
## GO:0015986 GO:0015986   ATP synthesis coupled proton transport     8/371
##            BgRatio       pvalue     p.adjust       qvalue
## GO:0042651 94/2476 1.332696e-08 1.959064e-06 1.753548e-06
## GO:0008137 19/2476 1.822145e-07 1.339277e-05 1.198780e-05
## GO:0019684 20/2476 4.505429e-07 2.207660e-05 1.976065e-05
## GO:0048038 21/2476 1.023581e-06 3.761661e-05 3.367044e-05
## GO:0042773 14/2476 2.908849e-06 8.552015e-05 7.654865e-05
## GO:0015986 10/2476 8.100531e-06 1.984630e-04 1.776432e-04
##                                                                                                                                                                                                                                                                                                     geneID
## GO:0042651 sll0519/sll0520/sll0522/sll0629/sll0689/sll1316/sll1322/sll1323/sll1324/sll1325/sll1326/sll1867/slr0228/slr0261/slr0782/slr0906/slr0927/slr1279/slr1280/slr1281/slr1291/slr1311/slr1329/slr1513/slr1623/slr1655/slr1796/sml0001/smr0006/smr0007/smr0008/smr0009/ssl2615/ssl3335/ssr1386/ssr3451
## GO:0008137                                                                                                                                                                                         sll0027/sll0520/sll0521/sll0522/sll1732/sll1733/slr0261/slr1279/slr1280/slr1281/slr1291/slr2007/slr2009
## GO:0019684                                                                                                                                                                                         sll0519/sll0520/sll0522/sll1867/slr0261/slr0906/slr0927/slr1279/slr1280/slr1281/slr1311/smr0006/ssr3451
## GO:0048038                                                                                                                                                                                         sll0519/sll0520/sll0521/sll0522/sll1732/slr0261/slr1279/slr1280/slr1281/slr1623/slr2007/slr2009/ssr1386
## GO:0042773                                                                                                                                                                                                                 sll0027/sll0522/sll1732/sll1733/slr1136/slr1291/slr2006/slr2007/slr2008/slr2009
## GO:0015986                                                                                                                                                                                                                                 sll1321/sll1322/sll1323/sll1324/sll1325/sll1326/slr1329/ssl2615
##            Count
## GO:0042651    36
## GO:0008137    13
## GO:0019684    13
## GO:0048038    13
## GO:0042773    10
## GO:0015986     8
go_functional_enrichment(CDS_result_dWT_dIF_0h, write=TRUE, path_up="output/enrichment/go_enrichment/dWT_dIF_0h_up.csv", path_down="output/enrichment/go_enrichment/dWT_dIF_0h_down.csv")
##                    ID Description GeneRatio  BgRatio       pvalue   p.adjust
## GO:0003677 GO:0003677 DNA binding      8/36 145/2423 0.0009785381 0.04892691
##                qvalue
## GO:0003677 0.04223165
##                                                                     geneID
## GO:0003677 sll0649/sll0856/sll1689/sll7106/slr0449/slr7049/slr7071/slr7092
##            Count
## GO:0003677     8
## [1] ID          Description GeneRatio   BgRatio     pvalue      p.adjust   
## [7] qvalue      geneID      Count      
## <0 rows> (or 0-length row.names)
go_functional_enrichment(CDS_result_dWT_dIF_1h, write=TRUE, path_up="output/enrichment/go_enrichment/dWT_dIF_1h_up.csv", path_down="output/enrichment/go_enrichment/dWT_dIF_1h_down.csv")
##                    ID                                      Description
## GO:0015979 GO:0015979                                   photosynthesis
## GO:0030089 GO:0030089                                    phycobilisome
## GO:0030096 GO:0030096 plasma membrane-derived thylakoid photosystem II
## GO:0042651 GO:0042651                               thylakoid membrane
## GO:0009523 GO:0009523                                   photosystem II
## GO:0051607 GO:0051607                                             <NA>
##            GeneRatio BgRatio       pvalue     p.adjust       qvalue
## GO:0015979    28/211 85/2473 6.415827e-11 8.148100e-09 7.766527e-09
## GO:0030089     9/211 24/2473 8.504556e-05 5.400393e-03 5.147495e-03
## GO:0030096     9/211 27/2473 2.427170e-04 7.707963e-03 7.347002e-03
## GO:0042651    19/211 94/2473 2.427705e-04 7.707963e-03 7.347002e-03
## GO:0009523     6/211 13/2473 3.703654e-04 9.407282e-03 8.966742e-03
## GO:0051607     5/211 10/2473 7.616605e-04 1.612181e-02 1.536683e-02
##                                                                                                                                                                                                                                     geneID
## GO:0015979 sll0047/sll0427/sll0629/sll1214/sll1398/sll1579/sll1580/slr0149/slr0506/slr1655/slr1739/slr1834/slr1835/slr1986/slr2051/slr2067/sml0001/sml0008/smr0004/smr0005/smr0006/smr0007/smr0008/ssl3093/ssr0330/ssr3304/ssr3383/ssr3451
## GO:0030089                                                                                                                                                         sll1579/sll1580/slr0149/slr1986/slr2051/slr2067/slr7023/ssl3093/ssr3383
## GO:0030096                                                                                                                                                         sll0427/sll1398/sll1867/slr0927/sml0001/smr0006/smr0007/smr0008/ssr3451
## GO:0042651                                                                         sll0047/sll0616/sll0629/sll1867/slr0927/slr1034/slr1655/slr1834/slr1835/slr1949/sml0001/sml0008/smr0004/smr0005/smr0006/smr0007/smr0008/ssl3335/ssr3451
## GO:0009523                                                                                                                                                                                 sll0047/sll1867/slr0927/smr0006/smr0008/ssr3451
## GO:0051607                                                                                                                                                                                         slr7016/slr7071/slr7092/ssr7017/ssr7093
##            Count
## GO:0015979    28
## GO:0030089     9
## GO:0030096     9
## GO:0042651    19
## GO:0009523     6
## GO:0051607     5
## [1] ID          Description GeneRatio   BgRatio     pvalue      p.adjust   
## [7] qvalue      geneID      Count      
## <0 rows> (or 0-length row.names)
kegg_functional_enrichment(CDS_result_dWT, write=TRUE, path_up="output/enrichment/kegg_enrichment/dWT_0h_1h_up.csv", path_down="output/enrichment/kegg_enrichment/dWT_0h_1h_down.csv")
## Reading KEGG annotation online:
## 
## Reading KEGG annotation online:
## [1] ID          Description GeneRatio   BgRatio     pvalue      p.adjust   
## [7] qvalue      geneID      Count      
## <0 rows> (or 0-length row.names)
##                ID               Description GeneRatio BgRatio       pvalue
## syn00190 syn00190 Oxidative phosphorylation    19/104  49/949 1.418619e-07
##              p.adjust       qvalue
## syn00190 6.525648e-06 6.421118e-06
##                                                                                                                                                           geneID
## syn00190 sll0223/sll0519/sll0520/sll0521/sll0522/sll1324/sll1325/sll1326/sll1327/sll1732/sll1733/slr0261/slr1279/slr1280/slr1281/slr1291/slr1330/slr1623/ssl2615
##          Count
## syn00190    19
kegg_functional_enrichment(CDS_result_dIF, write=TRUE, path_up="output/enrichment/kegg_enrichment/dIF_0h_1h_up.csv", path_down="output/enrichment/kegg_enrichment/dIF_0h_1h_down.csv")
## [1] ID          Description GeneRatio   BgRatio     pvalue      p.adjust   
## [7] qvalue      geneID      Count      
## <0 rows> (or 0-length row.names)
##                ID               Description GeneRatio BgRatio       pvalue
## syn00190 syn00190 Oxidative phosphorylation    23/172  49/949 1.760797e-06
## syn00195 syn00195            Photosynthesis    24/172  62/949 6.183843e-05
##              p.adjust       qvalue
## syn00190 9.860462e-05 9.638045e-05
## syn00195 1.731476e-03 1.692420e-03
##                                                                                                                                                                                                   geneID
## syn00190         sll0027/sll0519/sll0520/sll0521/sll0522/sll1322/sll1323/sll1324/sll1325/sll1326/sll1732/sll1733/slr0261/slr1136/slr1137/slr1279/slr1280/slr1281/slr1291/slr1329/slr1623/ssl2615/ssr1386
## syn00195 sll0427/sll0629/sll1316/sll1322/sll1323/sll1324/sll1325/sll1326/sll1398/sll1796/sll1867/slr0906/slr0927/slr1311/slr1329/slr1655/slr1828/sml0001/smr0006/smr0007/smr0008/ssl0020/ssl2615/ssr3451
##          Count
## syn00190    23
## syn00195    24
kegg_functional_enrichment(CDS_result_dWT_dIF_0h, write=TRUE, path_up="output/enrichment/kegg_enrichment/dWT_dIF_0h_up.csv", path_down="output/enrichment/kegg_enrichment/dWT_dIF_0h_down.csv")
## [1] ID          Description GeneRatio   BgRatio     pvalue      p.adjust   
## [7] qvalue      geneID      Count      
## <0 rows> (or 0-length row.names)
## --> No gene can be mapped....
## --> Expected input gene ID: sll1479,slr0301,slr1051,sll0745,slr0435,sll1213
## --> return NULL...
## NULL
kegg_functional_enrichment(CDS_result_dWT_dIF_1h, write=TRUE, path_up="output/enrichment/kegg_enrichment/dWT_dIF_1h_up.csv", path_down="output/enrichment/kegg_enrichment/dWT_dIF_1h_down.csv")
##                ID                       Description GeneRatio BgRatio
## syn00195 syn00195                    Photosynthesis     18/73  62/949
## syn00196 syn00196 Photosynthesis - antenna proteins      7/73  15/949
##                pvalue     p.adjust       qvalue
## syn00195 1.485540e-07 6.684929e-06 6.684929e-06
## syn00196 4.699591e-05 1.057408e-03 1.057408e-03
##                                                                                                                                                   geneID
## syn00195 sll0427/sll0629/sll1182/sll1398/sll1867/slr0927/slr1655/slr1739/slr1834/slr1835/sml0001/sml0008/smr0004/smr0005/smr0006/smr0007/smr0008/ssr3451
## syn00196                                                                                         sll1579/sll1580/slr1986/slr2051/slr2067/ssl3093/ssr3383
##          Count
## syn00195    18
## syn00196     7
## [1] ID          Description GeneRatio   BgRatio     pvalue      p.adjust   
## [7] qvalue      geneID      Count      
## <0 rows> (or 0-length row.names)

3.3 Gene Set Enrichment Analyses

dWT_0h_1h_go_gsea <- go_gsea(CDS_result_dWT, write=TRUE, path="output/enrichment/go_gsea/dWT_0h_1h_go_gsea.csv")
## preparing geneSet collections...
## GSEA analysis...
## Warning in fgseaMultilevel(...): For some pathways, in reality P-values are less
## than 1e-10. You can set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...
##                    ID                              Description setSize
## GO:0008137 GO:0008137 NADH dehydrogenase (ubiquinone) activity      19
## GO:0042773 GO:0042773 ATP synthesis coupled electron transport      14
## GO:0048038 GO:0048038                          quinone binding      21
## GO:0055114 GO:0055114              oxidation-reduction process     201
## GO:0019684 GO:0019684           photosynthesis, light reaction      20
## GO:0015977 GO:0015977                          carbon fixation      13
##            enrichmentScore       NES       pvalue     p.adjust      qvalues
## GO:0008137      -0.9218843 -2.682260 1.000000e-10 1.760000e-08 1.410526e-08
## GO:0042773      -0.9163969 -2.464829 2.748115e-09 2.418341e-07 1.938144e-07
## GO:0048038      -0.8397445 -2.465698 8.177258e-09 4.797325e-07 3.844746e-07
## GO:0055114      -0.4423196 -1.961674 3.030759e-07 1.333534e-05 1.068741e-05
## GO:0019684      -0.7784678 -2.261980 4.913748e-06 1.729639e-04 1.386194e-04
## GO:0015977      -0.8257245 -2.165586 1.894119e-05 5.556084e-04 4.452842e-04
##            rank                   leading_edge
## GO:0008137  120  tags=68%, list=2%, signal=67%
## GO:0042773   98  tags=64%, list=2%, signal=63%
## GO:0048038  120  tags=62%, list=2%, signal=61%
## GO:0055114  899 tags=28%, list=18%, signal=24%
## GO:0019684  120  tags=40%, list=2%, signal=39%
## GO:0015977  266  tags=62%, list=5%, signal=58%
dIF_0h_1h_go_gsea <- go_gsea(CDS_result_dIF, write=TRUE, path="output/enrichment/go_gsea/dIF_0h_1h_go_gsea.csv")
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
##                    ID                        Description setSize
## GO:0055114 GO:0055114        oxidation-reduction process     202
## GO:0042651 GO:0042651                 thylakoid membrane      94
## GO:0015979 GO:0015979                     photosynthesis      85
## GO:0006412 GO:0006412                        translation      89
## GO:0003735 GO:0003735 structural constituent of ribosome      55
## GO:0005840 GO:0005840                           ribosome      54
##            enrichmentScore       NES       pvalue     p.adjust      qvalues
## GO:0055114      -0.4276135 -2.112589 6.073483e-10 6.040886e-08 4.443953e-08
## GO:0042651      -0.5499966 -2.399096 6.864643e-10 6.040886e-08 4.443953e-08
## GO:0015979      -0.5494639 -2.369990 8.059436e-09 4.084434e-07 3.004697e-07
## GO:0006412      -0.5350417 -2.321256 9.282804e-09 4.084434e-07 3.004697e-07
## GO:0003735      -0.6065864 -2.385362 5.033700e-08 1.771863e-06 1.303463e-06
## GO:0005840      -0.6003029 -2.358980 1.347555e-07 3.952828e-06 2.907882e-06
##            rank                   leading_edge
## GO:0055114 1403 tags=42%, list=25%, signal=33%
## GO:0042651 1276 tags=59%, list=23%, signal=46%
## GO:0015979 1601 tags=61%, list=29%, signal=44%
## GO:0006412 1724 tags=71%, list=31%, signal=49%
## GO:0003735 1667 tags=80%, list=30%, signal=56%
## GO:0005840 1667 tags=80%, list=30%, signal=56%
dWT_dIF_0h_go_gsea <- go_gsea(CDS_result_dWT_dIF_0h, write=TRUE, path="output/enrichment/go_gsea/dWT_dIF_0h_go_gsea.csv")
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
##                    ID                                Description setSize
## GO:0004519 GO:0004519                      endonuclease activity      24
## GO:0006355 GO:0006355 regulation of transcription, DNA-templated     135
## GO:0003677 GO:0003677                                DNA binding     145
##            enrichmentScore      NES       pvalue    p.adjust     qvalues rank
## GO:0004519       0.7330815 2.018810 3.323554e-05 0.005782985 0.005527596  450
## GO:0006355       0.4386987 1.691613 4.164360e-04 0.028494167 0.027235804 1053
## GO:0003677       0.4361410 1.687431 4.912787e-04 0.028494167 0.027235804  739
##                              leading_edge
## GO:0004519 tags=46%, list=10%, signal=41%
## GO:0006355 tags=36%, list=24%, signal=28%
## GO:0003677 tags=27%, list=17%, signal=23%
dWT_dIF_1h_go_gsea <- go_gsea(CDS_result_dWT_dIF_1h, write=TRUE, path="output/enrichment/go_gsea/dWT_dIF_1h_go_gsea.csv")
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
##                    ID                        Description setSize
## GO:0015979 GO:0015979                     photosynthesis      85
## GO:0042651 GO:0042651                 thylakoid membrane      94
## GO:0006412 GO:0006412                        translation      89
## GO:0003735 GO:0003735 structural constituent of ribosome      55
## GO:0005840 GO:0005840                           ribosome      54
## GO:0030089 GO:0030089                      phycobilisome      24
##            enrichmentScore      NES       pvalue     p.adjust      qvalues rank
## GO:0015979       0.5607895 2.468703 3.516686e-10 6.189368e-08 5.441610e-08  864
## GO:0042651       0.4743434 2.098625 4.106142e-07 3.613405e-05 3.176857e-05 1375
## GO:0006412       0.4346562 1.920062 1.798293e-05 1.054999e-03 9.275406e-04 1634
## GO:0003735       0.4972601 1.983787 4.788636e-05 2.107000e-03 1.852446e-03 1563
## GO:0005840       0.4963100 1.980938 7.859357e-05 2.590323e-03 2.277378e-03 1563
## GO:0030089       0.6440629 2.125906 8.830648e-05 2.590323e-03 2.277378e-03  995
##                              leading_edge
## GO:0015979 tags=47%, list=16%, signal=40%
## GO:0042651 tags=46%, list=26%, signal=34%
## GO:0006412 tags=57%, list=31%, signal=40%
## GO:0003735 tags=65%, list=30%, signal=47%
## GO:0005840 tags=65%, list=30%, signal=46%
## GO:0030089 tags=58%, list=19%, signal=48%
dWT_0h_1h_kegg_gsea <- kegg_gsea(CDS_result_dWT, write=TRUE, path="output/enrichment/kegg_gsea/dWT_0h_1h_kegg_gsea.csv")
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
##                ID               Description setSize enrichmentScore      NES
## syn00190 syn00190 Oxidative phosphorylation      49      -0.7011069 -2.46859
##                pvalue     p.adjust      qvalues rank
## syn00190 8.115011e-09 4.950156e-07 4.014795e-07  896
##                            leading_edge
## syn00190 tags=59%, list=18%, signal=49%
dIF_0h_1h_kegg_gsea <- kegg_gsea(CDS_result_dIF, write=TRUE, path="output/enrichment/kegg_gsea/dIF_0h_1h_kegg_gsea.csv")
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
##                ID                           Description setSize enrichmentScore
## syn00190 syn00190             Oxidative phosphorylation      49      -0.6817810
## syn00195 syn00195                        Photosynthesis      62      -0.6048131
## syn03010 syn03010                              Ribosome      53      -0.6033970
## syn00240 syn00240                 Pyrimidine metabolism      25      -0.6715112
## syn01110 syn01110 Biosynthesis of secondary metabolites     285      -0.2879301
## syn03070 syn03070            Bacterial secretion system      14      -0.6774017
##                NES       pvalue     p.adjust      qvalues rank
## syn00190 -2.664246 2.962910e-10 1.807375e-08 1.341106e-08 1115
## syn00195 -2.495725 6.553581e-09 1.998842e-07 1.483179e-07 1608
## syn03010 -2.381502 8.731977e-08 1.775502e-06 1.317456e-06 1667
## syn00240 -2.213964 3.859505e-05 5.885745e-04 4.367335e-04  985
## syn01110 -1.483742 3.723699e-04 4.542913e-03 3.370928e-03 1417
## syn03070 -1.907194 2.445821e-03 2.486585e-02 1.845093e-02 1788
##                             leading_edge
## syn00190  tags=63%, list=20%, signal=51%
## syn00195  tags=73%, list=29%, signal=52%
## syn03010  tags=79%, list=30%, signal=56%
## syn00240  tags=52%, list=18%, signal=43%
## syn01110  tags=33%, list=26%, signal=26%
## syn03070 tags=100%, list=32%, signal=68%
dWT_dIF_0h_kegg_gsea <- kegg_gsea(CDS_result_dWT_dIF_0h, write=TRUE, path="output/enrichment/kegg_gsea/dWT_dIF_0h_kegg_gsea.csv")
## preparing geneSet collections...
## GSEA analysis...
## no term enriched under specific pvalueCutoff...
## [1] "nothing enriched"
dWT_dIF_1h_kegg_gsea <- kegg_gsea(CDS_result_dWT_dIF_1h, write=TRUE, path="output/enrichment/kegg_gsea/dWT_dIF_1h_kegg_gsea.csv")
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
##                ID                       Description setSize enrichmentScore
## syn00195 syn00195                    Photosynthesis      62       0.5836453
## syn03010 syn03010                          Ribosome      53       0.5039717
## syn00196 syn00196 Photosynthesis - antenna proteins      15       0.6968561
## syn00790 syn00790               Folate biosynthesis      17      -0.6603830
##                NES       pvalue     p.adjust      qvalues rank
## syn00195  2.367202 2.010784e-08 1.226578e-06 1.142972e-06 1044
## syn03010  1.971284 1.157646e-04 3.530821e-03 3.290152e-03 1563
## syn00196  2.037089 2.341332e-04 4.760709e-03 4.436208e-03  948
## syn00790 -1.835465 1.937164e-03 2.954175e-02 2.752812e-02  970
##                            leading_edge
## syn00195 tags=53%, list=20%, signal=43%
## syn03010 tags=68%, list=30%, signal=48%
## syn00196 tags=73%, list=18%, signal=60%
## syn00790 tags=47%, list=18%, signal=39%

3.4 Plots of GSEA

gseaplot2(dWT_0h_1h_go_gsea, geneSetID=1:6)

p <- gseaplot2(dWT_0h_1h_kegg_gsea, geneSetID =1)
p

ggsave("output/enrichment/Plots/dWT_0h_1h_KEGG.pdf", plot=p, width=15, height=12, units="cm")
p <- gseaplot2(dIF_0h_1h_kegg_gsea, geneSetID =1:6)
p

ggsave("output/enrichment/Plots/dIF_0h_1h_KEGG.pdf", plot=p, width=15, height=12, units="cm")
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn00190", 1) # https://www.kegg.jp/kegg-bin/show_pathway?syn00190/sll0026/slr2083/sll0223/sll1899/sll1262/sll1327/slr1622/sll1223/slr1279/slr1329/slr1137/sll0027/sll1326/ssr1386/sll0519/slr0261/sll1322/slr1136/slr1623/slr1280/sll1324/ssl2615/sll1325/sll1323/sll0522/sll0520/slr1281/sll0521/sll1733/slr1291/sll1732
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn00195", 2) # https://www.kegg.jp/kegg-bin/show_pathway?syn00195/sll1382/ssr0390/slr0343/sll0849/slr0342/sml0005/ssl2598/smr0005/slr0737/slr1645/slr1835/slr1643/smr0003/sll0819/sll1182/sll1327/slr1185/slr1834/smr0010/smr0004/sll1317/sll0427/smr0008/sll1796/slr1329/sll1316/slr1655/slr0906/sml0001/sll1326/slr1311/smr0006/smr0007/ssr3451/sll0629/slr1828/slr0927/sll1322/sll1324/ssl2615/sll1325/sll1323/sll1867/ssl0020/sll1398
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn03010", 3) # https://www.kegg.jp/kegg-bin/show_pathway?syn03010/slr1356/slr1678/sll1800/sll1260/sll1822/sll1801/sll1096/sll1821/sll1799/sll1808/ssr1604/sll1824/sll1816/sll1802/ssr1398/sll1804/sll1813/sll1806/ssl1784/sll1097/sll1807/ssr1399/ssl3437/sll1809/ssl3432/ssr0482/sll1746/sll1812/sll1817/ssl3436/ssr1736/sll1803/sll1805/sll1811/sll1244/sll0767/ssl1426/sll1743/sll1740/sll1745/slr1984/smr0011
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn00240", 4) # https://www.kegg.jp/kegg-bin/show_pathway?syn00240/sll1459/sll1108/slr1476/sll1635/slr1616/sll0144/slr1164/sll1258/sll0368/sll1035/sll1852/sll0370/slr0591
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn01110", 5) # https://www.kegg.jp/kegg-bin/show_pathway?syn01110/sll0902/sll0329/sll1558/slr1022/slr0738/slr0186/slr1289/slr1898/slr0049/slr0506/slr1887/slr1019/slr1808/slr0886/sll1444/slr1790/slr1867/slr0879/slr1055/sll1127/sll1468/sll1185/slr1226/sll1538/sll1172/sll1688/slr2136/slr0444/slr0293/slr1737/sll0513/sll1994/sll1899/sll0179/sll0901/sll1459/sll1108/sll0033/slr2072/sll1479/slr1349/slr0585/slr2130/sll1056/sll1747/slr0526/slr1124/sll1349/sll0158/sll1908/slr1176/sll1363/sll1091/slr0018/slr1510/slr1096/sll0807/sll0421/sll0480/slr1165/sll1342/slr0657/slr1743/slr0611/slr0985/slr1934/slr0966/slr0940/slr0116/slr1736/sll1815/slr1842/sll1797/slr1517/slr2094/slr0394/slr1945/sll0373/slr1843/sll0735/slr2023/sll0418/slr1030/sll0422/sll0018/sll1214/slr0749/slr0348/slr0012/slr0009/sll1852/sll0927/slr0493/sll0019/ssl2084
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn03070", 6) # https://www.kegg.jp/kegg-bin/show_pathway?syn03070/slr1471/sll1180/ssr3307/sll1814/sll1181/slr1531/sll0194/slr1046/sll0616/slr0774/slr0775/ssl2823/ssl3335


browseKEGGNew_3(dWT_0h_1h_kegg_gsea, "syn00190", 1) # https://www.kegg.jp/kegg-bin/show_pathway?syn00190/sll1899/slr2082/sll0026/sll1262/sll1322/slr1622/sll1484/slr1136/sll1323/slr1137/ssl2615/slr1330/slr1623/sll1327/sll1326/sll1324/sll1325/slr0261/slr1279/sll0223/sll0519/slr1280/sll0521/sll0520/sll0522/slr1281/slr1291/sll1733/sll1732
p <- gseaplot2(dWT_dIF_1h_kegg_gsea, geneSetID=1:4)
p

ggsave("output/enrichment/Plots/dWT_dIF_1h_KEGG.pdf", plot=p, width=15, height=12, units="cm")

p <- gseaplot2(dWT_dIF_1h_go_gsea, geneSetID=c(2,8,3))
p

ggsave("output/enrichment/Plots/dWT_dIF_1h_GO.pdf", plot=p, width=15, height=12, units="cm")

p <- gseaplot2(dWT_dIF_1h_kegg_gsea, geneSetID=c(1,3,2,4))
p

ggsave("output/enrichment/Plots/dWT_dIF_1h_KEGG.pdf", plot=p, width=15, height=12, units="cm")
browseKEGGNew_3(dWT_dIF_1h_kegg_gsea, "syn00195", 1) # https://www.kegg.jp/kegg-bin/show_pathway?syn00195/smr0005/slr1655/smr0004/sll0629/slr1834/slr1835/ssr3451/sll1182/smr0007/smr0006/smr0008/slr1739/sml0008/sll1867/sll0427/slr0927/sml0001/sll1398/sll1323/sll1382/slr1311/ssl0020/smr0010/sml0002/slr0737/sll0819/sll1322/slr0150/ssl0563/sll1325/slr0906/sll0849/sll1324
browseKEGGNew_3(dWT_dIF_1h_kegg_gsea, "syn03010", 2) # https://www.kegg.jp/kegg-bin/show_pathway?syn03010/sll1740/smr0011/slr1984/ssr0482/sll1811/sll1244/sll1743/ssl1784/ssr1604/sll1803/sll1812/ssl1426/sll1824/sll1813/sll1806/sll1807/ssr1736/ssl3437/sll1096/sll1805/ssl3432/sll1809/sll1810/ssl3436/sll1802/ssr1398/sll1804/slr1678/sll1799/ssr1399/sll1800/sll1808/sll1821/sll1097/sll1801/sll1260
browseKEGGNew_3(dWT_dIF_1h_kegg_gsea, "syn00196", 3) # https://www.kegg.jp/kegg-bin/show_pathway?syn00196/slr1986/slr2067/sll1579/sll1580/ssr3383/ssl3093/slr2051/sll1577/sll1578/slr0335/sll0928
browseKEGGNew_3(dWT_dIF_1h_kegg_gsea, "syn00790", 4) # https://www.kegg.jp/kegg-bin/show_pathway?syn00790/slr0527/slr1093/slr0078/slr0887/slr0900/slr0902/slr0901/slr0903

3.5 GSEA for RNA features

First, a data frame in which feature names are assigned to their feature type (CDS, ncRNA, …) has to be created and the respective info read in.

df_featureType <- as.data.frame(cbind(feature_type=as.character(features$type), feature_name=features$locus_tag))

feature_type_gsea <- function(DESEq2_dataframe, write=FALSE, path=""){
  locus_tags <- row.names(DESEq2_dataframe)
  geneList <- DESEq2_dataframe$log2FoldChange
  names(geneList) <- locus_tags
  geneList = sort(geneList, decreasing = TRUE)
  
  set.seed(42)
  go_gsea_object <- GSEA(geneList, TERM2GENE = df_featureType, seed=TRUE)
  print(head(go_gsea_object)[,1:10])
  
  if(write){
    write.csv(go_gsea_object, path)
  }
  return(go_gsea_object)
}

features_gsea_0h <- feature_type_gsea(CDS_result_dWT_dIF_0h, TRUE, "output/enrichment/GSEA_RNAfeatures_dWT-dIF_0h.csv")
## preparing geneSet collections...
## GSEA analysis...
## Warning in fgseaMultilevel(...): For some pathways, in reality P-values are less
## than 1e-10. You can set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...
##          ID Description setSize enrichmentScore       NES       pvalue
## 5UTR   5UTR        5UTR     272      -0.4118204 -2.119681 1.000000e-10
## tRNA   tRNA        tRNA      35       0.7360567  2.283511 5.012588e-07
## ncRNA ncRNA       ncRNA     223      -0.2909000 -1.443708 2.227298e-03
## 3UTR   3UTR        3UTR      25       0.6002223  1.745930 6.680201e-03
##           p.adjust qvalues rank                   leading_edge
## 5UTR  4.000000e-10      NA  846 tags=41%, list=20%, signal=35%
## tRNA  1.002518e-06      NA  652 tags=63%, list=15%, signal=54%
## ncRNA 2.969731e-03      NA  627 tags=28%, list=15%, signal=25%
## 3UTR  6.680201e-03      NA  363  tags=40%, list=8%, signal=37%
features_gsea_1h <- feature_type_gsea(CDS_result_dWT_dIF_1h, TRUE, "output/enrichment/GSEA_RNAfeatures_dWT-dIF_1h.csv")
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
##        ID Description setSize enrichmentScore      NES       pvalue
## tRNA tRNA        tRNA      38       0.6716029 2.459963 2.409659e-08
## 3UTR 3UTR        3UTR      29       0.5561780 1.910522 7.894904e-04
##          p.adjust      qvalues rank                   leading_edge
## tRNA 7.228978e-08 2.536484e-08  489  tags=53%, list=9%, signal=48%
## 3UTR 1.184236e-03 4.155212e-04 1062 tags=48%, list=20%, signal=39%
p <- gseaplot2(features_gsea_0h, geneSetID=1:4)
p

p <- gseaplot2(features_gsea_1h, geneSetID=1:2)
p

3.6 GSEA for base means

go_gsea_baseMeans <- function(DESEq2_dataframe, write=FALSE, path=""){
  locus_tags <- row.names(DESEq2_dataframe)
  geneList <- DESEq2_dataframe$baseMean
  names(geneList) <- locus_tags
  geneList = sort(geneList, decreasing = TRUE)
  
  set.seed(42)
  go_gsea_object <- GSEA(geneList, TERM2GENE = term_to_gene, TERM2NAME=term_to_name, seed=TRUE)
  tryCatch({
    print(head(go_gsea_object)[,1:10])
  
  if(write){
    write.csv(go_gsea_object, path)
  }
  return(go_gsea_object)}, error=function(e){
    print("nothing enriched")
  })
}

kegg_gsea_baseMean <- function(DESEq2_dataframe, write=FALSE, path=""){
  locus_tags <- row.names(DESEq2_dataframe)
  geneList <- DESEq2_dataframe$baseMean
  names(geneList) <- locus_tags
  geneList = sort(geneList, decreasing = TRUE)
  
  set.seed(42)
  kegg_gsea_object <- gseKEGG(geneList, organism="syn", minGSSize=10, pvalueCutoff = 0.05, seed=TRUE)
  tryCatch({
    print(head(kegg_gsea_object)[,1:10])
  
  if(write){
    write.csv(kegg_gsea_object, path)
  }
  return(kegg_gsea_object)}, error=function(e){
    print("nothing enriched")
  })
}

3.7 Without taking width into account

When width of features is not taken into account, there is no specific enrichment of certain terms in more highly expressed features. One exception is the GO term “transposase activity”, which is depleted when analyzing the comparison of dWT and dIF at time point 0h.

go_gsea_baseMeans(CDS_result_dIF)
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".
## Warning in fgseaMultilevel(...): There were 2 pathways for which P-values were
## not calculated properly due to unbalanced (positive and negative) gene-level
## statistic values. For such pathways pval, padj, NES, log2err are set to NA. You
## can try to increase the value of the argument nPermSimple (for example set it
## nPermSimple = 10000)
## no term enriched under specific pvalueCutoff...
## [1] "nothing enriched"
go_gsea_baseMeans(CDS_result_dWT)
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There were 2 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple =
## 10000)
## no term enriched under specific pvalueCutoff...
## [1] "nothing enriched"
go_gsea_baseMeans(CDS_result_dWT_dIF_0h)
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".
## leading edge analysis...
## done...
##                    ID          Description setSize enrichmentScore       NES
## GO:0004803 GO:0004803 transposase activity      11      -0.7504063 -2.409325
##                  pvalue    p.adjust     qvalues rank
## GO:0004803 1.316033e-05 0.002289897 0.002271888 1087
##                               leading_edge
## GO:0004803 tags=100%, list=25%, signal=75%
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism     UNKNOWN 
## #...@setType      UNKNOWN 
## #...@geneList     Named num [1:4318] 209545 165602 81721 65565 59371 ...
##  - attr(*, "names")= chr [1:4318] "ncl1680" "ncr0700" "ncr0075" "sll1578" ...
## #...nPerm     
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...1 enriched terms found
## 'data.frame':    1 obs. of  11 variables:
##  $ ID             : chr "GO:0004803"
##  $ Description    : chr "transposase activity"
##  $ setSize        : int 11
##  $ enrichmentScore: num -0.75
##  $ NES            : num -2.41
##  $ pvalue         : num 1.32e-05
##  $ p.adjust       : num 0.00229
##  $ qvalues        : num 0.00227
##  $ rank           : num 1087
##  $ leading_edge   : chr "tags=100%, list=25%, signal=75%"
##  $ core_enrichment: chr "sll1861/sll1710/ssl3649/sll1716/slr0800/slr2096/slr1903/slr5040/sll1792/sll1157"
## #...Citation
##   Guangchuang Yu, Li-Gen Wang, Yanyan Han and Qing-Yu He.
##   clusterProfiler: an R package for comparing biological themes among
##   gene clusters. OMICS: A Journal of Integrative Biology
##   2012, 16(5):284-287
go_gsea_baseMeans(CDS_result_dWT_dIF_1h)
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There were 2 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple =
## 10000)
## no term enriched under specific pvalueCutoff...
## [1] "nothing enriched"
kegg_gsea_baseMean(CDS_result_dIF)
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".
## no term enriched under specific pvalueCutoff...
## [1] "nothing enriched"
kegg_gsea_baseMean(CDS_result_dWT)
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".
## no term enriched under specific pvalueCutoff...
## [1] "nothing enriched"
kegg_gsea_baseMean(CDS_result_dWT_dIF_0h)
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".
## no term enriched under specific pvalueCutoff...
## [1] "nothing enriched"
kegg_gsea_baseMean(CDS_result_dWT_dIF_1h)
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".
## no term enriched under specific pvalueCutoff...
## [1] "nothing enriched"

3.8 When taking width of features into account

norm_counts <- counts(ddsMat_CDS, normalized=TRUE)
mean_norm_counts <- apply(norm_counts, 1, mean)
mean_norm_counts_CDS <- subset(mean_norm_counts, names(mean_norm_counts) %in% CDS_features$locus_tag)
CDS_feat_df <- as.data.frame(CDS_features)
row.names(CDS_feat_df) <- CDS_feat_df$locus_tag
mean_norm_counts_CDS_width <- mean_norm_counts_CDS/CDS_feat_df[names(mean_norm_counts_CDS),]$width
mean_norm_counts_CDS_width <- sort(mean_norm_counts_CDS_width, decreasing=TRUE)
set.seed(42)
go_gsea_baseMeans_width <- GSEA(mean_norm_counts_CDS_width, TERM2GENE=term_to_gene, TERM2NAME = term_to_name, nPermSimple=10000, seed=TRUE)
## preparing geneSet collections...
## GSEA analysis...
## Warning in fgseaMultilevel(...): There were 3 pathways for which P-values were
## not calculated properly due to unbalanced (positive and negative) gene-level
## statistic values. For such pathways pval, padj, NES, log2err are set to NA. You
## can try to increase the value of the argument nPermSimple (for example set it
## nPermSimple = 100000)
## leading edge analysis...
## done...
set.seed(42)
kegg_gsea_baseMeans_width <- gseKEGG(mean_norm_counts_CDS_width, organism="syn", minGSSize=10, pvalueCutoff = 0.05, nPermSimple=10000, seed=TRUE)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
go_gsea_baseMeans_width[,1:10]
##                    ID                 Description setSize enrichmentScore
## GO:0015979 GO:0015979              photosynthesis      86       0.9254062
## GO:0042651 GO:0042651          thylakoid membrane      96       0.8945113
## GO:0018298 GO:0018298 protein-chromophore linkage      31       0.9451642
## GO:0008272 GO:0008272           sulfate transport      12      -0.5812176
##                  NES       pvalue     p.adjust      qvalues rank
## GO:0015979  1.565034 4.810606e-07 8.466667e-05 8.102074e-05  189
## GO:0042651  1.509690 1.417049e-04 1.247003e-02 1.193305e-02  196
## GO:0018298  1.621643 8.657711e-04 4.930474e-02 4.718157e-02  127
## GO:0008272 -1.811787 1.120562e-03 4.930474e-02 4.718157e-02 1547
##                               leading_edge
## GO:0015979   tags=56%, list=5%, signal=54%
## GO:0042651   tags=50%, list=5%, signal=49%
## GO:0018298   tags=45%, list=3%, signal=44%
## GO:0008272 tags=100%, list=42%, signal=58%
write.csv(go_gsea_baseMeans_width, "output/enrichment/GO_GSEA_baseMeans_perWidth.csv")
kegg_gsea_baseMeans_width[,1:10]
##                ID    Description setSize enrichmentScore      NES       pvalue
## syn00195 syn00195 Photosynthesis      63       0.9294204 1.580491 2.778544e-05
##             p.adjust     qvalues rank                  leading_edge
## syn00195 0.001722697 0.001722697  195 tags=70%, list=5%, signal=67%
write.csv(kegg_gsea_baseMeans_width, "output/enrichment/KEGG_GSEA_baseMeans_perWidth.csv")
p <- gseaplot2(go_gsea_baseMeans_width, geneSetID=c(1,2,3,4))
p

ggsave("output/enrichment/Plots/GO_GSEA_baseMeans_perWidth.pdf", plot=p, width=15, height=12, units="cm")
p <- gseaplot2(kegg_gsea_baseMeans_width, geneSetID=c(1))
p

ggsave("output/enrichment/Plots/KEGG_GSEA_baseMeans_perWidth.pdf", plot=p, width=15, height=12, units="cm")

Session info

## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 11 (bullseye)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8    
##  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=de_DE.UTF-8   
##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] forcats_0.5.1               stringr_1.4.0              
##  [3] dplyr_1.0.5                 purrr_0.3.4                
##  [5] readr_1.4.0                 tidyr_1.1.3                
##  [7] tibble_3.1.0                tidyverse_1.3.0            
##  [9] rtracklayer_1.50.0          enrichplot_1.10.2          
## [11] clusterProfiler_3.18.1      RColorBrewer_1.1-2         
## [13] pheatmap_1.0.12             ggpubr_0.4.0               
## [15] ggrepel_0.9.1               ggplot2_3.3.3              
## [17] DESeq2_1.30.1               SummarizedExperiment_1.20.0
## [19] Biobase_2.50.0              MatrixGenerics_1.2.1       
## [21] matrixStats_0.58.0          GenomicRanges_1.42.0       
## [23] GenomeInfoDb_1.26.2         IRanges_2.24.1             
## [25] S4Vectors_0.28.1            BiocGenerics_0.36.0        
## [27] knitr_1.31                 
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1             shadowtext_0.0.7         backports_1.2.1         
##   [4] fastmatch_1.1-0          plyr_1.8.6               igraph_1.2.6            
##   [7] splines_4.0.4            BiocParallel_1.24.1      digest_0.6.27           
##  [10] htmltools_0.5.1.1        GOSemSim_2.16.1          viridis_0.5.1           
##  [13] GO.db_3.12.1             fansi_0.4.2              magrittr_2.0.1          
##  [16] memoise_2.0.0            openxlsx_4.2.3           Biostrings_2.58.0       
##  [19] annotate_1.68.0          graphlayouts_0.7.1       modelr_0.1.8            
##  [22] colorspace_2.0-0         rvest_1.0.0              blob_1.2.1              
##  [25] haven_2.3.1              xfun_0.22                crayon_1.4.1            
##  [28] RCurl_1.98-1.2           jsonlite_1.7.2           scatterpie_0.1.5        
##  [31] genefilter_1.72.1        survival_3.2-7           glue_1.4.2              
##  [34] polyclip_1.10-0          gtable_0.3.0             zlibbioc_1.36.0         
##  [37] XVector_0.30.0           DelayedArray_0.16.2      car_3.0-10              
##  [40] abind_1.4-5              scales_1.1.1             DOSE_3.16.0             
##  [43] DBI_1.1.1                rstatix_0.7.0            Rcpp_1.0.6              
##  [46] viridisLite_0.3.0        xtable_1.8-4             foreign_0.8-81          
##  [49] bit_4.0.4                httr_1.4.2               fgsea_1.16.0            
##  [52] ellipsis_0.3.1           pkgconfig_2.0.3          XML_3.99-0.5            
##  [55] farver_2.1.0             sass_0.3.1               dbplyr_2.1.0            
##  [58] locfit_1.5-9.4           utf8_1.1.4               labeling_0.4.2          
##  [61] tidyselect_1.1.0         rlang_0.4.10             reshape2_1.4.4          
##  [64] AnnotationDbi_1.52.0     munsell_0.5.0            cellranger_1.1.0        
##  [67] tools_4.0.4              cachem_1.0.4             cli_2.3.1               
##  [70] downloader_0.4           generics_0.1.0           RSQLite_2.2.3           
##  [73] broom_0.7.5              evaluate_0.14            fastmap_1.1.0           
##  [76] yaml_2.2.1               bit64_4.0.5              fs_1.5.0                
##  [79] tidygraph_1.2.0          zip_2.1.1                ggraph_2.0.5            
##  [82] xml2_1.3.2               DO.db_2.9                rstudioapi_0.13         
##  [85] compiler_4.0.4           curl_4.3                 ggsignif_0.6.1          
##  [88] reprex_1.0.0             tweenr_1.0.1             geneplotter_1.68.0      
##  [91] bslib_0.2.4              stringi_1.5.3            highr_0.8               
##  [94] lattice_0.20-41          Matrix_1.3-2             vctrs_0.3.6             
##  [97] pillar_1.5.1             lifecycle_1.0.0          BiocManager_1.30.10     
## [100] jquerylib_0.1.3          data.table_1.14.0        cowplot_1.1.1           
## [103] bitops_1.0-6             qvalue_2.22.0            R6_2.5.0                
## [106] gridExtra_2.3            rio_0.5.26               MASS_7.3-53.1           
## [109] assertthat_0.2.1         withr_2.4.1              GenomicAlignments_1.26.0
## [112] Rsamtools_2.6.0          GenomeInfoDbData_1.2.4   hms_1.0.0               
## [115] grid_4.0.4               rmarkdown_2.7            rvcheck_0.1.8           
## [118] carData_3.0-4            ggforce_0.3.3            lubridate_1.7.10